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A self-consistent version of the Thermal Random Phase Approximation 
(TSCRPA) is developed within the Matsubara Green's Function (GF) formal- 
ism. The TSCRPA is applied to the many level pairing model. The normal 
phase of the system is considered. The TSCRPA results are compared with 
the exact ones calculated for the Grand Canonical Ensemble. Advantages of 
the TSCRPA over the Thermal Mean Field Approximation (TMFA) and the 
standard Thermal Random Phase Approximation (TRPA) are demonstrated. 
Results for correlation functions, excitation energies, single particle level den- 
sities, etc., as a function of temperature are presented. 

I. INTRODUCTION 

Pairing properties of finite Fermi systems such as ultrasmall metallic grains have recently 
received a great deal of attention. This has been spurred by a series of spectacular experi- 
ments of Ralph, Black and Tinkham [1]. In order to correctly describe pairing properties it 
has been recognized that the finiteness of the systems (grains) needs to consider quantum 
fluctuations, good particle number, number parity, etc. seriously, since the coherence length 
may be of the order of the system size. The situation for metallic grains has in the mean- 
while been well described in several review articles [2,3] (see also [4]). Another system where 
the finiteness is at the forefront of the theoretical investigation since several decades is the 
superfluid atomic nucleus. As a matter of fact many of the theoretical tools such as particle 
number projection, even-odd effects, number parity, blocking effect, particle - particle Ran- 
dom Phase Approximation (pp-RPA), etc. have first been developed in nuclear physics [5] 
before finding their application to finite systems of condensed matter. Also the schematic 
pairing model with which we will mostly deal in this paper, namely the Picket Fence Model 
(PFM), whose exact solution has been found by Richardson and Sherman [6], has essentially 
been developed in the context of nuclear physics for the description of deformed superfluid 
nuclei. For finite condensed matter systems an early theoretical description was proposed by 
Muhlschlegel, Scalapino, and Denton [7] using the Static Path Approximation (SPA) to the 
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partition function. This work stayed rather singular for a long time but the SPA has recently 
been applied successfully to the PFM both in the condensed matter [8] and nuclear [9,10] 
contexts. A further standard method to treat quantum fluctuations namely the well known 
RPA has quite extensively been used for nuclear systems [5,11] but equally for condensed 
matter problems [12]. 

In this work we will further elaborate on the RPA approach. We indeed have recently 
had quite remarkable success with a self consistent extension of the pp-RPA, which we called 
Self-Consistent RPA (SCRPA), by reproducing very accurately groundstate and excitation 
energies of the PFM [13] at zero temperature. This formalism was also developed indepen- 
dently by Ropke and collaborators who called it Cluster - Hartree - Fock (CHF) [14]. Such 
type of generalization of the RPA theory grew out of the works of K.-J. Hara and D. Rowe [15] 
several decades ago. Shortly afterwards the theory was rederived using the method of many 
body Green functions [16]. The success of the theory motivates us to develop the SCRPA 
formalism also for the finite temperature case and to study the thermodynamic properties of 
the BCS Hamiltonian using the PFM as an example. For the extension of SCRPA to finite 
temperature we use the Matsubara Green functions approach [17]. It appeared that the ap- 
proximation scheme is very effective in treating two-body correlations in the particle-particle 
(pp) channel as well as the Pauli principle effects. We should mention that we will work with 
real particles and not with quasiparticles what should limit our approach to temperatures 
above the critical temperature T c (i.e. to the normal phase). However, as we will see below, 
the definition of T c in SCRPA is not so clear and we will be able to continue our calculation 
quite deeply into the superfluid regime. 

We organize the paper in the following way. In Section 2, the approach is outlined in 
general. Then, in Section 3, the formalism is applied to the PFM. A comparison with the 
exact solutions as well as with the results of other approximations is made in Section 4. 
Section 5 is devoted to comparison with other recent works. In section 6, we will summarize 
the results and draw some conclusions. In an appendix a variant for the calculation of the 
occupation numbers is proposed. 

II. GENERAL FORMALISM 

In treating a finite many-body system at finite temperature, it is convenient to use the 
grand canonical ensemble although it violates the number conservation. With the definition 

K = H - iiN 

the grand partition function and statistical operator read 

Z G = = Tr(e^ K ) 

where (3 = 1/T. Then for any Schrodinger operator A a the modified Heisenberg picture can 
be introduced 
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A a (r) = e Kr A a e- Kr 

and the temperature (or the Matsubara) Green's Function (GF) is defined as [17] 
Gtf = - (T T A a (T)A+(T>)) = -Tr [e^ K ^T T e^ A a e^ K A^' K 



-Tr 



,T T e TK A a e-^ K A + a e- T ' K 



(1) 



Here, the brackets ( ) mean the thermodynamic average; T T is a r ordering operator, which 
arranges operators with the earliest r (the closest to —f3) to the right. 
Let us consider the two-body Hamiltonian 



12 4 1234 



(2) 



where a, a + are fermion annihilation and creation operators; t± 2 and V1234 = V1234 — ^1243 

are the kinetic energy and the antisymmetrized matrix element of the two-body interaction. 
The Green's function G T ~^ for an arbitrary operator A+ obeys the following equation of 



motion: 



d 



A,, At]) - (T T [A a ,K] T AUt')) = 



In this expression it is possible to split the effective Hamiltonian Ti^ into an instantaneous 
and a dynamic (frequency dependent) part [18] 



a/3 



= E S r-r' ( [[A a , K] , A+\ )-{T T [A a , K] T [K, A+ 

0' L 

- n ap d T-T> + n a/3 



In the approximation of the instantaneous effective Hamiltonian i.e. neglecting T , the 

Dyson equation for the two-body Matsubara GF can be written as 



dr 



(3) 



In the treatment of two particle correlations let us specify the arbitrary operator A a as 
Ak ± k 2 — afciOfc 2 - I 11 ^is case the Dyson equation (3) takes the following form in the frequency 
representation 

r SCRPA (A\ 
2^ rL k 1 k 2Pl p 2 (jr p 1 p 2 k'k' > V*) 



%LO n^k x k 2 k'^ - iV /ci/c 2 fcm 



P\P2 



where, in supposing that the single particle density matrix is diagonal in the basis used (this 
is for example the case inhomogeneous matter): 



ak 1 ak2, a t 2l a t ll ) — <W 2 fcifc 2 (l - n k 2 , - n^,) 



(5) 
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and 

n kLk[k' 2 = E ( [K°fe , K] , «] ) Np^k'^ ( 6 ) 



PlP2 



Here Sk 1 k 2 k' 1 k' 2 is the antisymmetrized Kronecker symbol and nk = (a^ak/ are the single 
particle occupation numbers which can be found from the single-particle Matsubara GF 

G\ k , = (T T o fc (r)ot(0)> 

as 

n k = (ata k ) = lim G T k (7) 

T — >0 

In general the single-particle Matsubara Green's function G T k ~ T ' obeys the following Dyson 
equation: 

+ £*) Glk' = S(t) + J dnMl^G? (8) 
or in the frequency representation 

Gt = Gl + G° k M k Gt n (9) 

where 

Gl = (10) 

lU n — Ek 

Here Ek contains already the usual (instantaneous) mean field so that denotes only the 
dynamical part of the mass operator. 

Now the problem is to find an approximation for the mass operator consistent with 
the SCRPA. A solution to this problem has been proposed in [18], which goes via the two 
body T - matrix representation of the single particle mass operator [19], evaluating T within 
SCRPA. 

Let us add at this point a word of physical interpretation of the mean-field operator 
(6) [13]. A quick look allows to realize that it contains no higher than two body correla- 
tion functions and therefore for their determination, with (4), one obtains a selfconsistency 
problem. Furthermore one can consider the nucleus as a gas of zero point pair fluctuations. 
These fluctuations create their own mean field, i.e. one pair fluctuation moves in the average 
potential created by all the other pair fluctuations. This average pair fluctuation field is 
graphically represented in Fig. 1. It gives rise, as usual, to a nonlinear problem. Of course, 
the single particle mean field introduced in (9) and further developed below is coupled to 
the selfconsistent pair potential. This is the deeper meaning of of (6). 




FIG. 1. Schematic representation of the first order self-energy of a pair fluctuation. Exchange 
terms are not presented. The full black dot is the interaction. 
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III. APPLICATION TO THE PICKET FENCE MODEL 



The model consists of an equidistant multilevel pairing Hamiltonian with each level two 
fold degenerate, i.e. only spin up/down fermions of one kind can occupy one level. The 
corresponding Hamiltonian is given by 

H = J2 e kN k -GJ2 p t p k (11) 

fc=l i,k=l 

with 

N k = 4c k + c^c^, 

K = 44 (12) 

where k means the time reversed of k, single particle energies are e k = ke — A, with level 
spacing e chosen to be equal to 1, and f2 stands for the number of levels. The chemical 
potential A will be chosen such as to conserve the average number of particles N = f2 of the 
system. The operators defined in (12) form an SU(2) algebra for each level j and obey the 
following commutation relations 

'Pi,*?] =Mi-^), 

[P j ,N k \=25j k P j , (13) 



A. SCRPA equations 

To study the model at finite temperature we define in analogy to (1) the following set of 
two-body Matsubara GFs 

C% = -(r T P,-(r)pt(0)), 

where 

1 

3 



\/< u-^-i > 

Applying the instantaneous approximation for the mass operator we obtain the expressions 
for the two body SCRPA GF's: 

^u Jn Gf RPA = 5 3^ + Y.^G^ p ^ (14) 

k 

with 

«(»>-2* (r G V rP*P ^ ^ < (1 - N,)(l - N t ) > 

- 1 [ e > + 71^1^ E < p > p >' > ) - G v /|< 1 - J v J >< 1 -iv t >| (15) 

To find the correlation functions of the form < (1 — Nj)(l — N k ) > we will use the 
following approximation: 
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when j ^ k it is a simple factorization procedure, which has turned out to be accurate in 
the zero temperature limit: 



< (i _ Nj )(i _ Nk ) >=< 1 — Nj >< 1 — N k > 
but when j = k we use the following exact relation 

< (1 - Nj) 2 >=< 1 - Nj > +2 < P+Pj > 



(16a); 



(166) 



which can easily be obtained taking into account that n 2 = rij and rijUj = P^Pj (here 
rij = CjCj). It should be noted that in [13] the factorization (16a) was also used for the 
diagonal part (16b) and quite accurate results were obtained. We will show below that 
with (16b) one obtains still improved results. As shown in [13] it is possible to avoid above 
approximation. However, this is at the cost of a considerable numerical complication. We 
refrain from this here because it brings only very little improvement of results. 
With this ansatz a particle-particle RPA-like equation is obtained 



SCRPA 



Gj\DjDi\ 



z-C, (z-CMz-G) 



x 



1 + GS r% 

k Z — Ofc 



;i7) 



where 



Z = lUJ r , 



and 



D, =< 1 - N > 



G 



C j = 2[e j -Gn j + — J2<P+P j ,> 



3 j'ft 



From this one easily can find the excitation spectrum of the model in equating the denomi- 
nator of (17) to zero 

1 + GEA=0 (18) 

k Z — Ofc 

Knowing the poles of the Green's function (17), one can write down its spectral representation 
(we here give it as a function of imaginary time), with the corresponding residua: 



G T P1P2 = (r) JD P1 D P2 X»X»e- E ^ (1 + n B (E,)) + Y»Y£e E ^n B {E,) 



+ 



+6(-t) ^D Pl D P2 
-Gl lh2 = 9 (r) JD^D 



X^ 2 e- E ^n B {E,) + Y£Y£e E " r (1 + n B {E,)) 
Y&Yfer E >" (1 + n B {E,)) + X^X^n B {E, 



+ 



(19a) 



+9 (-r) ^D hl D h2 [Y/;X 2 e- E ^n B (E,) + X^Xfc^ (1 + n B {E,)) 



(196) 
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-G T ph = -G\ v = 9 (r) ^\D p D h \ [x^e^ (1 + n B {E,)) + Y»X^n B {E») 



+ 



+9 (-r) ^\D p D h \ [i;y;e^n B (£,) + Y^e^ (1 + n B {E,)) 

1 



(19c) 



where the index p refers to the states above Fermi level and the index h to the ones below. 
The following amplitudes were introduced in these formulas: 



\C P \ — E^ 



F = 



\Ch\ + E^ 



with 



\Ch\ — E^ 

i+gj: 



GD 

F yfj- — Y p 



y^-h — 



F~ 2 = — 

" dz 



\C P \ + E^ ^ 



k z — C k 

These amplitudes obey the usual normalization conditions 



(20a) 
(206) 
(20c) 



z=E„ 



E*W + EW = *„ 

V h 



(21) 



h V 

Two-body correlation functions can be obtained from the Green's function (19) as follows 



< P+P» >= -G^- = ^D pi D p2 £ X^ 2 n B (E,) + £ Y£Y£ (1 + n B (E,)) 



^D hl D h2 Y£Y h >B(E») + E (! + n B{E,)) (22) 



< P p + P h >= ~GW°~ = yl\D P D h \ (Y X ^ Y hME,) + E Y?K U + ME,)) 



B. Occupation numbers in the SCRPA 

In order to close the set of the SCRPA equations, it is necessary to find the so far 
unknown occupation numbers n k =< c\c k > . For this, we should find the single particle 
Green's function G T k consistent with the SCRPA scheme. As discussed in sect. II, the single 
particle mass operator M k has in general the exact representation in terms of the two body 
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T-matrix [19] and then an appropriate approximation for the G\~ can be obtained. It consists 
in using the mass operator M k calculated through the T-matrix found in the framework of 
SCRPA. As the relation between the T-matrix and the sum of the all irreducible Feynman 
graphs in the pp-chanel is also known then the following expression for the single particle 
mass operator can be obtained [18]: 

M k = Gj2 Gf'^Gt^H^ (23) 

k\k2 

Hfak i s expressed through the effective Hamiltonian (15) without the disconnected part 

ntl = nfl - 2S kkl e k (24) 

where e k is defined below in (30). 

In addition to this transparent scheme there also exists an additional consistency require- 
ment [20]. It follows from the possibility to calculate the average value of the Hamiltonian 
< H > in two ways. On the one hand one has the following relation between the single 
particle Green's function (9) and < H > [17]: 

< H >= ~\ , >!"V E \l - A K" T,) + G r'') < 25 > 

k 



On the other hand there exists the straightforward calculation of < H > through the two 
body Green's functions (22) 

< H >= ]Te fc < N k > -G £ < P+P k2 > 

k k\k2 

= - T + 2 E e p<^>- 2 EE ^/GDpF, {n B {E»)X» + (1 + n B (E,))Y») (26) 
v v v 

where we used the particle-hole symmetry of the model [11] and reduced sums over p and h 
only to the one over the particle states. 

The additional consistency condition lies in the requirement that both expressions (25) 
and (26) should give exactly the same results. This only is satisfied if the single particle GF 
is expanded to first order in the renormalized single particle mass operator (23) (one may 
verify that this is in analogy to the standard RPA scheme, i.e. the standard RPA average 
energy is obtained via (25) using a single particle GF with only perturbative renormalization 
from the RPA-modes): 

Gp = G° p + G° p M p SCRPA G° p (27) 

with 

M p CRPA = ^M P) (28) 



D° — 1— f — f- 

L Jp J] 



p 



f> = ttW- < 29) 
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and 



£ = e -Gf 

I Op U p 



(30) 



Finally we get the following expression for the SCRPA single particle mass operator 



M SCRPA. 



M: 



SCRPA 



Jgd p 

£)0 2^ r V 



p n 



S$(f P + n B (E,)) f Sj$(l-f p + n B (E„)) 



iw 

iuj n + e p — E^ 



+ 



iuj n + e p + E^ 



(31) 



where 



sS = -EWp' + E^?V, 

p' h' 

p' /i' 

The corresponding single particle occupation numbers n p , found from (7), is the following 



Tip <C C p Cp ^> fp ~\~ 



GD r 



E^ 



£)0 

P M 



r^gj^-g K(^) + l)D° + / p 2 
p (2e p -^) + p (2e p + E,) 



GD r 



fv (1 - f P ) E ^ *pU + n B {E,) + -f p + n B (E,) 



D° 

V V 



(32) 



Let us demonstrate now that using the single particle Green's function (27) with the 
mass operator (31) and occupation numbers (32) in the calculation of the average energy 
(25) indeed leads to the equation (26). At first one finds the derivative of the single particle 
GF (27) 



£pfp~^~ 



GD r , 



D° 



Emu-/*) 



c(D ( £ p - E^riBjE^ (2) (e p + E^)(l + n B (E^) 

>->,,„ — — — r o. m 



MP 



(2% - EJ* 



'HP 



(2e p + EJ* 



~^~£pfp 



q (i) fp + n B {E^) (2) 1 - fp + nsjEn 

MP (Or _ P ^2 + (Or 



+/fe J ,/ p (l " /„) 



{2e p -Etf ^ (2e p + £, 

ff(2) l-/ p + 7l B (^) 



2e v + E u 



/p 



» 2e p -E, + » 2 £p + ^ JJ 



-£ p n p -\- ^ ' 



GDpF^ 



qnaiEjiq - H s(2) n B (E,)D° p + f t 



2e p — E^ 



vp 



2e p + E^, 



Inserting this expression in (25) we obtain 

< H >= lim V 

2 t'-t^o+ V 
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= —r + 2 E( e p + £ pK - 2 E y/GDpF* {n B {E,)X» + (1 + n B (^))y/) + G E ^/p = 

= — r + 4 S e ^ - 2 S E \[gd^ {"BiEjx* + (i + 

P P M 

This is exactly equal to (26). 

The system of the SCRPA equations is fully closed now. Together with (17), (19) and 
(22) this represents a self-consistent problem for pair fluctuations. 

We want to indicate at this point that the above way to determine the single particle 
occupancies is not the only possibility. In the Appendix we will give another variant which, 
however, yields results close to the ones with the method of this section. The non uniqueness 
of the occupation numbers reflects the fact that with the truncated ansatz (3), at zero 
temperature, no corresponding ground state wave function can be found, as explained in 
[13]. For a wave function to exist, the ansatz (3) must be extended. It can, however, be 
shown that the correction terms are small [21]. 



C. Exact statistical treatment of the PFM 



For an exact statistical treatment of the Picket Fence Model we have to find all exact 
eigenvalues and eigenstates of the Hamiltonian (11). Since singly-occupied levels do not 
participate in the pair scattering, eigenstates can be classified according to the number of 
unpaired particles S (seniority). There are Cq_ s = srmz^y different multiplets of this type, 
each of dimension C ^lf and degeneracy 2 s . If we define the following set of basis states for 

2 

each multiplet: 

\{s i ,N i }> (33) 

where s$ = iVj = 1 for singly-occupied levels and s$ = 0, iVj = or 2 for remaining lev- 
els Si = S), the Hamiltonian matrix will have the following diagonal and off-diagonal 
elements 

({s l ,N i }\H\{s i ,N t }) = ^2(e k -\)+ E 

fees feen-s 

( s h N jn s i fe 2, Sj fi, s jn N jn \H\ s^N^, s jk 0, s it 2, s jn N jn ) = -G (35) 

The exact eigenvalues and eigenstates can be calculated by diagonalization of this matrix 
in each multiplet. The exact grand canonical average < A > of any operator can then be 
obtained with the help of the grand partition function Zq = Tre~^ H and the statistical 
operator p G = Z G 1 e~ /3H as 

< A >= Tr[Ap G ] (36) 



(e fc - A) - T (4 - MO 



N, 



(34) 
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IV. RESULTS AND DISCUSSION 



In order to check the accuracy of our theory and of the different approximations schemes 
we first calculate the average energy of the system < H > as a function of the particle 
number N and temperature T. The results of the different calculations are presented in 
Figures 2-4. Calculations were made for a value of the pairing constant G which is smaller 
but close to the critical value G cr at T = 0. The phase transition from the normal to the 
superfluid phase occurs in the system when G — > G cr . We compare the SCRPA results with 
the exact ones for the Grand Canonical Ensemble (GCE) as well as results of the standard 
Thermal RPA (TRPA) and Thermal MFA (TMFA). One can see that when the number of 
levels Q (and number of particles N) increases the description of the intrinsic energy becomes 
better and at Q = N = 10 the TSCRPA results practically coincide with the exact ones. 
Especially the last case will be considered below more carefully. 

Let us now come to the discussion of the behavior of the excitation energies. The depen- 
dence of the excitation energies of the addition mode (see [13]) as a function of G is shown 
in Fig. 5 at zero temperature. The SCRPA (solid lines) is compared with the standard RPA 
calculations (dashed lines) and the exact ones (open circles). Increasing the interaction con- 
stant, the lowest energy in the RPA goes to zero and at G cr ~ 0.33 the collapse takes place 
which is connected with the transition from the normal to the superfluid phase. At finite 
temperature (see Fig. 6 where the dependence of the lowest addition mode is presented as 
a function of G at T = 0.5) this collapse occurs at a higher value of the interaction constant 
(G cr — 0.43) what is due to the reduced intensity of the residual interaction because of the 
thermal factors. This collapse is absent in the exact calculations at zero temperature and 
also in the SCRPA calculations at zero and finite temperatures. It is also remarkable that 
the SCRPA yields a rise of all excitation energies with increasing G in contrast with RPA and 
in very good agreement with the exact results. This comes from the fact that in the PFM 
with the Kramer's degeneracy of levels the Pauli repulsion is extremely strong overruling the 
original attractive interaction. In this model, therefore, standard RPA gives qualitatively 
wrong result. 

We next consider the behavior of the system near the phase transition point. To make 
distinctions between different results more apparent we not only show the full intrinsic energy 
< H > but also the correlation energy E corr which is defined as 

E corr =<H>-<H> (37) 

where < H > is the average energy calculated in Mean Field Approximation. In Figures 
7 and 8, the average energy < H > and correlation energy E corr as a function of T are 
displayed for the interaction constant G = 0.4 (at T = this value of G is larger than 
G cr ~ 0.33). With increasing T the mean field rearrangement occurs and the system goes 
from the superfluid phase to the normal one at T c ~ 0.38. Note, that within the TRPA 
the lowest excitation energy uji vanishes when T — > T c , whereas within the TSCRPA uo\ 
stays finite. For both the correlation energy and the intrinsic energy the TSCRPA gives 
more precise results as compared to the other approximations. It is remarkable that the 
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TSCRPA results are accurate down to practical zero temperature, in spite of the fact that 
within standard BCS theory one enters the superfluid regime. A quasiparticle formulation 
of SCRPA [22] will only be necessary for stronger G values driving the system more deeply 
into the symmetry broken phase. 

To analyze the region near the phase transition point in more detail, the heat capacity 
is calculated as a partial derivative of the intrinsic energy with respect to T 

C = ^ (38) 

The results are shown in Fig. 9. The TRPA and TMFA give discontinuities of C v at T c (we 
recall again that our results are obtained using a normal fluid approach and not transforming 
to quasiparticles). The heat capacity calculated in the TSCRPA has some kink near T c but 
has no discontinuities and is quite similar to the exact result through out the whole range 
of temperature. 

Nevertheless, the TSCRPA and also the exact solution feel the phase transition to the 
superfluid phase. It can be seen in Fig. 8 where both the TSCRPA and the exact correlation 
energies show a depression near T = 0.38. This originates from strong pair fluctuations 
leading to the BCS phase transition in TMFA with the critical temperature T RCS = 0.38 for 
G = 0.4. However, one notices (see also [23]) that the sharp phase transition of mean field 
is in reality completely smeared out and only a faint, though clearly visible, signal survives. 

To investigate the formation of such fluctuations in more detail, it is useful to consider 
the spectral function [17] 

1 



A(k,u) 



<^k 



~ u k 



(39) 



which includes two-body correlations through the self-consistent treatment of the mass op- 
erator Mj? CRPA (uj n ). Based on the spectral function the density of states can be calculated 
as [24] 

N(u) — ^2,A(k,u) (40), 

k 

The results of the calculations of N(u) with G^ of (27) in (39) are shown in Fig. 10 at 
different values of T and for G = 0.4. It is clearly seen that the distance between the 
two quasiparticle peaks around the Fermi energy ef (oj = 0) increases with decreasing 
temperature. This process sets in even above the BCS transition temperature T c s cs = 0.38. 
This rarefaction of the level density around e F above T c is not avoid of similarity with the 
situation in high T c - superconductors where a so-called 'pseudo gap' in the level density 
appears already above T c [25]. This 'pseudo gap' also is often attributed to a decrease in 
the level density around ef due to pair fluctuations [24,26]. Apparently it is a quite generic 
feature that pair correlations diminish the density of levels around Bp whereas particle-hole 
correlations give rise to an increase. 

In order to make the temperature dependence of the gap more transparent let us introduce 
an effective (or canonical) gap which recently was proposed in Eq. (22) of Ref. [3] 



A = GJY, (<^> - (4ck) (etc,)) (41) 

V ik 
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In the BCS approximation, the effective gap A coincides with the usual grand canonical 
BCS gap 

A = Gj2(Pk)Bcs 

k 

The dependencies of the effective gap on the interaction constant G at zero temperature and 
on temperature T at G = 0.4 are shown in Figures 11 and 12. The SCRPA results give a 
very good description of the gap at zero and non zero temperatures. It is clearly seen that 
the SCRPA and exact calculations do not display the phase transition at the point where 
the BSC gap vanishes. Notice that in Fig. 11 the SCRPA result deteriorates for values of G 
well beyond the critical value. In this regime a quasiparticle generalization of the SCRPA is 
necessary [22] because one enters deeply in the superfluid region. 

V. COMPARISON WITH OTHER WORKS 

The PFM has recently widely been used for the study of quantum pair fluctuations at 
finite temperature both in the context of nuclear physics [9,10] and of ultrasmall metallic 
grains [8,4]. In both fields the SPA approach where one additionally takes into account 
number parity projection and quantum (RPA) fluctuation around mean field was employed. 

To compare our results with the above mentioned formalism let us introduce the rel- 
evant energy scales. These are the average level spacing 5 and the BCS energy A = 
Q8/[2 sinh(5/G)] [2,3]. The properties of the system described by the pairing Hamiltonian 
can be calculated as universal functions of the single scaling parameter 8/ A. As long as the 
grain is not too small 8 « A, the fluctuation region around T c is narrow, and the mean 
field (BCS) description of superconductivity is appropriate. It is the case in ref. [10], where 
the PFM is investigated with characteristic values of 8/ A ~ 0.2. The result of that work 
shows that the domain where the parity projection is important lies in a small region near 
T c . At temperatures higher than T c the role of the fluctuation is decreased and usual SPA 
becomes rather good in reproducing the exact canonical results. 

When the size of the system is decreased, fluctuations start to smear the normal - super- 
conducting transition. The finite level spacing suppresses the BCS gap and when 8 becomes 
of the order of A, the fluctuation region becomes of order T c and the BCS description of 
superconductivity breaks down even at zero temperature. Returning to our calculations, 
we can see that SCRPA yields the best results for 8/ A > 1 (see Fig. 11, 12). This region 
corresponds to ultrasmall grains where strong pairing fluctuations are dominant. In this 
sense our results can be compared (at least qualitatively) with the results of [8]. To make 
this comparison more accurate we calculated average energy and specific heat for a system 
with 50 fermions and interaction constant G = 0.127 and G = 0.255, what corresponds to 
6/ A = 50 and 6/ A = 1 respectively. In general our results correlate well with [8]. From Fig. 
13, where C v is displayed as a function of T/8, we can see that when 5/ A — > oo the specific 
heat approaches a linear dependence with temperature while when 8 ~ A a bump structure 
arises at low temperature which is a sign of the presence of strong pairing fluctuations. As 
it has been seen from our previous discussion, the SCRPA gives a better description with 
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increasing particles number. And while we did not perform exact GCE calculations for the 
system with 50 particles, we can expect from our studies above and in [13] that our result 
should be very close to the exact one. 

In conclusion of this section we can say that the results of TSCRPA are at least as good 
as the ones of SPA with extensions. However, contrary to SPA, TSCRPA has no problem at 
low temperature and excitation energies and correlation functions can be calculated directly 
as a function of temperature. 

VI. CONCLUSION 

In this work we generalized our recent work [13] on the multilevel pairing model (PFM) 
within the SCRPA approach to finite temperature (TSCRPA). The PFM has been recognised 
to account in many respects for the physics of (superconducting) metallic nano-grains. In our 
context SCRPA can be viewed as a self-consistent mean field theory for pair fluctuatuations. 
The results at T = in [13] are in very close agreement with the exact ones obtained from 
the Richardson procedure [5]. It is therefore an important issue to also exploit SCRPA at 
finite T and to assess its accuracy with respect to exact results in this case. Our comparison 
is mostly done for the case of ten levels with 10 particles where it is still of some ease to 
establish the exact partition function. We, however, also considered with TSCRPA the case 
of 50 particles in 50 levels assuming that the results are of equal quality or even better as 
the ones obtained for 10 particles. We base our studies on the Matsubara 1 and 2 particle 
Green's functions which allows us to calculate correlation and excitation energies, specific 
heat, level densities, etc. It can be considered as a general advantage of our approach that 
all these quantities are directly accessible in the whole range of temperatures and coupling 
constants. For the latter this holds in this work only true for interaction values not driving 
the system deeply into the superfluid regime, since in this work we only have been working 
with normal particles and not with quasiparticles. For G » G cr we have to employ the 
Self Consistent Quasiparticle RPA (SCQRPA). It has recently been demonstrated in the two 
level pairing model that also SCQRPA gives very promising results [22]. The quality of our 
results for the above mentioned quantities turns out to be excellent and it does not fail in any 
qualitative nor quantitative aspect. Most of the time the agreement with the exact solution 
is within the couple of percent level. One particularly interesting feature of our investigation 
is the fact that we achieved to calculate the single particle Green's function consistently 
within TSCRPA. This enabled us to give, for the first time, for the PFM the evolution of the 
single particle level density with temperature. The construction of the exact solution for this 
quantity is very cumbersome and we refrained from working this out here. However, backed 
with the positive experience for all other above mentioned quantities we believe that also the 
level density is reasonably accurate. The interesting aspect of our calculation is that with 
decreasing temperature the density of single particle states around the Fermi level decreases 
even above the critical temperature as defined by BCS - theory. It is suggestive to see this 
feature in analogy to the appearance of a so called pseudogap in high T c superconductors 
where also a depression in the level density is observed approaching T c from above [24,25]. 
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It would be interesting to attempt an experimental verification with metallic nanograins of 
our prediction that indeed the density of levels rarefies with decreasing temperature already 
in the non superfluid regime. 

We also gave a short comparison of TSCRPA with results at finite T obtained with other 
approaches like the static path approximation (SPA) to the partition function. Though 
the results seem generally comparable, we think that TSCRPA is more versatile, giving 
direct access to correlation functions, level densities, excitation energies, etc. in the whole 
temperature range, quantities which are otherwise difficult to obtain. 

In the present work we restricted ourselves to values of the coupling which are below 
or slightly above the critical value. In the future we shall elaborate on the SCRPA for 
quasiparticles at finite T (TSCQRPA) which will allow us to consider the system deeply in 
the superfluid phase and to study the transition from one phase to the other in more detail. 
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APPENDIX: OCCUPATION NUMBER - VARIANT 



Let us give another possible way to find the occupation numbers at finite temperature. 
Before to complete this task we firstly derive the SCRPA equations with the Equation of 
Motion method and find some useful relations between phonon amplitudes at zero tempera- 
ture. To attain this let us introduce the pair addition and removal operators (phonons) [13] 

as 

p h 

Qi = J2 x kP h -T,Yp x Pp m 

h p 

and apply the variational procedure where all expectation values are found with respect to 
the vacuum of phonons (A1,A2). 



< 



5Qu, ft, Qt 



>= OJv < 



5Qu,Qt 



> 



(A3) 



If we use the factorisation procedure (16) the following system of equations for the phonon 
amplitudes is obtained 



p' 
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where 



~~ 25 pp / [e p ] + 2G5 pp / 



^ D 



- n n 



— G\J DpD p r — 25p P 'Ep + A 



pp' 



Bph — —G\j—D v D h 



Chh' = 2(W[eft - G] + 2G5, 



hh> 



l^h D h 



n h 



Bh P — —G\j—D h Dp 
where the single-particle energies are introduced as 



— G\J D h D h i = — (25hh>£h + Chh') 

{A5) 



E p — e p , Eh — —eh + G 
Thereafter the expressions for the phonon amplitudes are obtained as 

rfl _ - Y, P > X^App, + J2h' Yh'Bph' 



p> h> 



XL 



2e p - 



(A6) 



{Ala) 



Y£{2e h + E,) = -Y: X$B W + £ Y&C-, 



hh> 



Y jx = - Hp' X^Bhpi + Eh' Yh' c hh' ^ A?b ^ 



2E h + 



P ' h> 

Due to the particle-hole symmetry, the removal mode satisfy exactly the same equations. It 
means that 

yf _ Y x =^ — V A= ^ ( AR\ 

^■p — ^h=N-p+li I h — 1 p=N-h+l 

Below we will use the notation ji for both modes. 

Returning to the single particle occupation numbers n k let us remind that in the picket 
fence model the following exact relation between one body operator N k and two body oper- 
ator P k Pk is verified for any non singly-occupied level k: 



N k = 2P+P k 



{A9) 



Taking the average of both parts of this relation with respect to the SCRPA vacuum state 
we obtain 

< N p >= 2D p J2{Y p n 2 0410) 
On the other hand, using relations (A7) we can rewrite this expression as 

1 



<Np>=2D p J2 



h p' 



h p' 



PP 



{All) 



r(E, + 2e p ) 

It is easy to show that the fraction in this expression can be expressed through the GFs as 
1 



{Ep + 2e p f 



- (( lim + J J dt^dhGf-^O^ - i 1 ) e - i ^+^)(*i-*i) G o( t ' 1 - t ') (A12) 
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where 

Gf = -i< T tCp (t)cJ-(ti) >= -i8(t - tJe-^-V (A13) 
If we now take into account the spectral representation of the pair operator P p 



P p (t) = D p 



(A14) 



and define the antichronological (time reversed) single particle GF as 

Gf- tl] = -i< Ttc+ftcpfa) >= -i9(t! - t^'V-V (Alb) 

we can pass to the expression for the sp GF which gives occupation numbers which are 
consistent within the frame of the SCRPA with exact relation (A9) 

= Of'*') + J J dUdt^Gf-^M^^Gf'^ (A16) 
where M p SCRPA is a single particle mass operator 

Mp CRPA = £ ^"^G^ (A17) 

fcifc2 

and Ti^l is the renormalised effective Hamiltonian which has the following form in terms of 
RPA matrixes A pp r and B ph 

H { pp l = Jd~ p A pp ,-^= (A18a) 



n;> = p p B ph ^= (Am) 



517(0) 

,/;/(:,, 

From (A16) and (A17) we then can define the occupation numbers as usual. In general, the 
occupation numbers from (32) and (A16), (A17) have slightly different values. This is due to 
the fact that the ansatz (Al), (A2) for the RPA operators is too restricted for a groundstate 
fulfilling <5|0 >= to exist. It can, however, be shown [21] that the necessary corrections 
to (Al), (A2) are small. Anyhow, the differences in occupation numbers obtained from the 
two methods are a measure of the importance of the terms neglected in (Al), (A2). The 
difference of the mass operators (23) and (A17) is that in the latter both vertices are dressed, 
whereas in (23) one vertex remains at the unrenormalized value (G). 

To find occupation numbers < n p > at finite temperature we adopt the expressions 
obtained for single particle GF at zero temperature (eq.(Al6)-(Al8)). To do it, it is necessary 
to change all zero temperature GFs to the Matsubara ones and use for the vertices the 
renormalised effective Hamiltonian (15) 



H% = JD p n^^== (A19) 



where 

7(0) _ o,(0) 
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and 



£fc — <2fc — Gf k , 
fk = 



1 + e £ ^ 

Then we get for the single particle Matsubara GF 



Mt-t') _ nO(T-r') 



T— Tl) 



? 0(r(-n)^7(0) ^(tJ-tiW ) 



\ - ^0(r(-ri)7r7(0) ~( 
2^ H-pk^k 



nO(ri-T') 



where 



Gf = - e -^(— ') [0 (r - r') (1 - /„) - 61 (r' - r) /J , 



(420) 



Taking this integral in the limit r — r' — ► we obtain the occupation numbers n p in the 
SCRPA at finite temperature 

n p =< c; Cp >= U + D P Y: (Y») 2 [n B (E,) (1 - 2f p ) + (1 - f p ) 



(X;) 2 [n B (E,) (1 - 2f p ) - {fp) 



-fp(l-f P )PD p 



J2(Y»)\l + n B (E,)-f p ) (2e p + E,) 



-f p (l-f p )PD p 



]T(x;) 2 (n B (2e p -E„) 



(A21) 



One should note here one thing about the direct use of the exact relation (A9) for the 
definition of the < N p >. In reality the identity (A9) is only valid for the collective subspace 
(spanned by the non singly-occupied levels) or, in other words, for levels which have partial 
seniority = 0. But when we work in the grand canonical ensemble, seniority of the level 
s k is not a good quantum number, since averaging procedure over GCE mixes all seniorities. 
This fact is reflected in the above expression (A21) where we can see that level occupation 
numbers < N p >= n p + n p due to thermal factors f p (Fermi-Dirac distribution) are not equal 
to < 2P+P p >. 

Before coming to a numerical example let us make a further comment. Knowing the 
occupation numbers as a function of the amplitudes X and Y as in (A21) or (32), a natural 
idea would be to use this to express the ground state energy entirely as a function of X 
and Y. Then minimizing under the constraint of the normalization conditions (21) leads 
to an equation determining X, Y amplitudes. Such a procedure has been proposed in the 
past by Jolos and Rybarska [27]. However, already in our example where the occupation 
numbers are exactly known (at zero temperature) as a function of X, Y amplitudes, we 
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have checked that this leads to worse results than with the approach advanced in the main 
text. In other examples the functional n p [X, Y] may only be known approximately and then 
the minimization of the ground state leads to deteriorated results, as we also have checked. 
In the Green's function approach where one first calculates excitation energies, i.e. energy 
differences, before calculating the ground state energy, such uncertainties like the precise 
knowledge of n p [X,Y] are minimized and the solution of the SCRPA equations is therefore 
much more stable. 

Let us now compare the results based on the two different definitions of the occupation 
numbers (eq. (32) and (A21)). We will denote results with eq. (32) by TSCRPA and results 
obtained with the second definition (eq. (A21)) by TSCRPAl. All calculations are made 
for the system with 10 particles on 10 levels. Results for the correlation energy (37) as a 
function of the coupling constant G are given in Tab. 1 and Tab. 2 for T = and T = 1, 
respectively. In these tables we denote results obtained with (A21) and (26) by TSCRPAl-1 
and the ones obtained with (A21) and (25) by TSCRPA1-2. We can see that TSCRPA 
systematically gives underestimated results with respect to the exact ones. Using definition 
(A21) we then can see that TSCRPAl-1 at zero temperature slightly overestimates the exact 
ground state energy (this corresponds to the results given in [13]) while TSCRPAl-2 always 
underestimates it. At finite temperature all approximations (TSCRPA and TSCRPAl) give 
close results which underestimate the exact ones. In Tab. 3 and 4 we show excitation energies 
of the first addition mode for T = and T = 1. In both cases, TSCRPA and TSCRPAl, 
excitation energies grow with increasing G for zero and nonzero temperatures reproducing 
quite well the exact results. We can see that for the calculated quantities (correlation and 
excitation energies) TSCRPAl-1 is slightly closer to the exact results than TSCRPA. 
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Table 1. Correlation energy as function of G at zero temperature. 



G 


Exact 


TSCRPA1-1 


TSCRPA1-2 


TSCRPA 


0.10 


-0.036 


-0.037 


-0.036 


-0.036 


0.20 


-0.167 


-0.169 


-0.153 


-0.159 


0.30 


-0.435 


-0.445 


-0.342 


-0.379 


0.33 


-0.551 


-0.564 


-0.408 


-0.461 


0.34 


-0.593 


-0.608 


-0.430 


-0.489 


0.35 


-0.638 


-0.654 


-0.453 


-0.518 


0.36 


-0.685 


-0.702 


-0.476 


-0.548 


0.40 


-0.896 


-0.917 


-0.569 


-0.670 



Table 2. Correlation energy as function of G at T — 1. 



G 


Exact 


TSCRPA1-1 


TSCRPA1-2 


TSCRPA 


0.10 


-0.030 


-0.029 


-0.018 


-0.029 


0.20 


-0.142 


-0.126 


-0.082 


-0.122 


0.30 


-0.372 


-0.299 


-0.208 


-0.295 


0.33 


-0.476 


-0.367 


-0.259 


-0.365 


0.34 


-0.515 


-0.392 


-0.278 


-0.391 


0.35 


-0.556 


-0.418 


-0.297 


-0.417 


0.36 


-0.600 


-0.444 


-0.317 


-0.445 


0.40 


-0.803 


-0.559 


-0.403 


-0.569 



Table 3. Excitation energy of the first addition mode as function of G at T = 0. 



G 


Exact 


TSCRPA1 


TSCRPA 


0.10 


1.001 


1.001 


1.001 


0.20 


1.005 


1.012 


1.012 


0.30 


1.014 


1.049 


1.053 


0.33 


1.018 


1.068 


1.074 


0.34 


1.019 


1.075 


1.081 


0.35 


1.022 


1.082 


1.089 


0.36 


1.023 


1.089 


1.098 


0.40 


1.031 


1.123 


1.136 
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Table 4. Excitation energy of the first addition mode as function of G at T = 1.0. 



G 


TSCRPA1 


TSCRPA 


0.10 


1.017 


1.030 


0.20 


1.075 


1.126 


0.30 


1.175 


1.281 


0.40 


1.299 


1.476 


0.41 


1.312 


1.497 


0.42 


1.324 


1.518 


0.43 


1.337 


1.539 


0.44 


1.349 


1.560 


0.45 


1.362 


1.581 
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FIG. 2. The average energy < H > as a function of the temperature for VL = N = 2 and G = 0.9. 
The exact results - open circles; the TMFA results - dotted line; the TRPA results - dashed line 
and the TSCRPA results - solid line. 
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FIG. 3. The average energy < H > as a function of the temperature for f2 = N = 4 and G = 0.5. 
For notation, see Fig.2 



23 



0.5 1 1.5 2 

T 



FIG. 4. The average energy < H > as a function of the temperature for 0, = N = 10 and 
G = 0.33. For notation, see Fig.2. 
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FIG. 5. Excitation energies as a function of the interaction constant G for f2 = 10 and T = 0. 
Notations: the exact results - open circles; TRPA results - dashed lines and the TSCRPA results 
- solid lines. 
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FIG. 6. The first excitation energy E\ as a function of the interaction constant G for VL 
and T = 0.5. Notations: TRPA results - dashed line, the TSCRPA results - solid line. 
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FIG. 7. The average energy < H > as a function of the temperature for 17 = N = 10 and 
G = 0.4. For notation, see Fig.2. 
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FIG. 8. The correlation energy E corr as a function of the temperature for f2 
G = 0.4. For notation, see Fig.2. 
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FIG. 10. The density of states as a function of the frequency u for Q = N = 10, G = 0.4 and 
T = 0.05, T = 0.25, T = 0.45 and T = 1.05. 
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FIG. 11. The effective gap A as a function of the interaction constant G for 17 = N = 10 and 
T = 0. The exact results - open circles; the BCS results - dotted line; the RPA results - dashed 
line and the SCRPA results - solid line. 
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FIG. 12. The effective gap A as a function of the temperature T calculated for £1 = N = 10 and 
G = 0.4. For notation, see Fig. 11. 
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FIG. 13. The specific heat C v as a function of the temperature calculated in the thermal SCRPA 
for Q = N = 50. Solid line corresponds to G = 0.128 (S/A = 50) and dotted line corresponds to 
G = 0.256 (S/A = 1). 
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